Factor 5
MOFAobject = load_model('../data/MOFA_models/MOFA2_noScale_train2020_21.hdf5')
MOFAobject
## Trained MOFA with the following characteristics:
## Number of views: 6
## Views names: geneExp ab cytokineL cell_freq tcellPol tcellAct
## Number of features (per view): 6650 31 14 39 1 2
## Number of groups: 1
## Groups names: group1
## Number of samples (per group): 158
## Number of factors: 15
plot_variance_explained(MOFAobject, x="view", y="factor")

plot_factor(MOFAobject,
factor = 5,
color_by = "IgG_PT"
)

plot_weights(MOFAobject,
#view = "geneExp",
factor = 5,
nfeatures = 20, # Number of features to highlight
scale = T, # Scale weights from -1 to 1
abs = F # Take the absolute value?
)

plot_weights(MOFAobject,
view = "geneExp",
factor = 5,
nfeatures = 10, # Number of features to highlight
scale = T, # Scale weights from -1 to 1
abs = F # Take the absolute value?
)

plot_weights(MOFAobject,
view = "cell_freq",
factor = 5,
nfeatures = 10, # Number of features to highlight
scale = T, # Scale weights from -1 to 1
abs = F # Take the absolute value?
)

plot_weights(MOFAobject,
view = "ab",
factor = 5,
nfeatures = 10, # Number of features to highlight
scale = T, # Scale weights from -1 to 1
abs = F # Take the absolute value?
)

plot_weights(MOFAobject,
view = "cytokineL",
factor = 5,
nfeatures = 10, # Number of features to highlight
scale = T, # Scale weights from -1 to 1
abs = F # Take the absolute value?
)

plot_top_weights(MOFAobject,
#view = "geneExp",
factor = 5,
nfeatures = 40
)

plot_top_weights(MOFAobject,
view = "cell_freq",
factor = 5,
nfeatures = 40
)

plot_data_heatmap(MOFAobject,
view = "geneExp", # view of interest
factor = 5, # factor of interest
features = 60, # number of features to plot (they are selected by weight)
# extra arguments that are passed to the `pheatmap` function
cluster_rows = TRUE, cluster_cols = FALSE,
show_rownames = TRUE, show_colnames = FALSE
)

plot_data_scatter(MOFAobject,
view = "geneExp", # view of interest
factor = 5, # factor of interest
features = 11, # number of features to plot (they are selected by weight)
add_lm = TRUE, # add linear regression
color_by = "Meta.infancy_vac"
)

plot_data_scatter(MOFAobject,
view = "cytokineL", # view of interest
factor = 5, # factor of interest
features = 11, # number of features to plot (they are selected by weight)
add_lm = TRUE, # add linear regression
color_by = "Meta.infancy_vac"
)

factors <- get_factors(MOFAobject, factors = "all")
lapply(factors,dim)
## $group1
## [1] 158 15
weights <- get_weights(MOFAobject, views = "all", factors = "Factor5")
a= weights$geneExp
ggplot(a, aes(x=Factor5)) +
geom_density()

a['Expr.ENSG00000277632.1', ]
## [1] 0.2673311
a = a[order(a[,1]), ]
head(a)
## Expr.ENSG00000121807.5 Expr.ENSG00000127951.6 Expr.ENSG00000091106.18
## -0.7597985 -0.6950563 -0.6691434
## Expr.ENSG00000165389.6 Expr.ENSG00000138061.11 Expr.ENSG00000150337.13
## -0.6136788 -0.6036396 -0.5989425
tail(a)
## Expr.ENSG00000182809.10 Expr.ENSG00000161835.10 Expr.ENSG00000105655.18
## 0.4987166 0.5021038 0.5158566
## Expr.ENSG00000177595.17 Expr.ENSG00000118503.14 Expr.ENSG00000178038.16
## 0.5259299 0.5611030 0.5694787
tail(names(a), n=200)
## [1] "Expr.ENSG00000172366.19" "Expr.ENSG00000007376.7"
## [3] "Expr.ENSG00000128159.11" "Expr.ENSG00000134285.10"
## [5] "Expr.ENSG00000175115.11" "Expr.ENSG00000167685.14"
## [7] "Expr.ENSG00000165915.13" "Expr.ENSG00000148362.10"
## [9] "Expr.ENSG00000181026.14" "Expr.ENSG00000169592.14"
## [11] "Expr.ENSG00000169429.10" "Expr.ENSG00000141582.14"
## [13] "Expr.ENSG00000197728.9" "Expr.ENSG00000064545.14"
## [15] "Expr.ENSG00000061273.17" "Expr.ENSG00000139626.15"
## [17] "Expr.ENSG00000204348.9" "Expr.ENSG00000010322.15"
## [19] "Expr.ENSG00000141232.4" "Expr.ENSG00000111674.8"
## [21] "Expr.ENSG00000161267.11" "Expr.ENSG00000066117.14"
## [23] "Expr.ENSG00000196155.12" "Expr.ENSG00000160285.14"
## [25] "Expr.ENSG00000132002.7" "Expr.ENSG00000136802.11"
## [27] "Expr.ENSG00000132382.14" "Expr.ENSG00000071894.16"
## [29] "Expr.ENSG00000137070.17" "Expr.ENSG00000106628.10"
## [31] "Expr.ENSG00000102225.15" "Expr.ENSG00000177606.6"
## [33] "Expr.ENSG00000125912.10" "Expr.ENSG00000114779.19"
## [35] "Expr.ENSG00000102901.12" "Expr.ENSG00000115268.9"
## [37] "Expr.ENSG00000112149.9" "Expr.ENSG00000167526.13"
## [39] "Expr.ENSG00000261236.7" "Expr.ENSG00000068831.18"
## [41] "Expr.ENSG00000025770.18" "Expr.ENSG00000142765.17"
## [43] "Expr.ENSG00000118046.14" "Expr.ENSG00000171130.17"
## [45] "Expr.ENSG00000204351.11" "Expr.ENSG00000149743.13"
## [47] "Expr.ENSG00000167615.16" "Expr.ENSG00000136573.12"
## [49] "Expr.ENSG00000110987.8" "Expr.ENSG00000087074.7"
## [51] "Expr.ENSG00000111678.10" "Expr.ENSG00000186654.20"
## [53] "Expr.ENSG00000179532.12" "Expr.ENSG00000198520.10"
## [55] "Expr.ENSG00000125652.7" "Expr.ENSG00000187189.10"
## [57] "Expr.ENSG00000162496.8" "Expr.ENSG00000108963.17"
## [59] "Expr.ENSG00000165271.16" "Expr.ENSG00000100003.17"
## [61] "Expr.ENSG00000165802.22" "Expr.ENSG00000167747.14"
## [63] "Expr.ENSG00000076928.17" "Expr.ENSG00000099821.13"
## [65] "Expr.ENSG00000171425.9" "Expr.ENSG00000109906.13"
## [67] "Expr.ENSG00000115687.13" "Expr.ENSG00000152518.7"
## [69] "Expr.ENSG00000171222.10" "Expr.ENSG00000112511.17"
## [71] "Expr.ENSG00000198932.12" "Expr.ENSG00000064547.13"
## [73] "Expr.ENSG00000074370.17" "Expr.ENSG00000185101.12"
## [75] "Expr.ENSG00000124181.14" "Expr.ENSG00000142102.15"
## [77] "Expr.ENSG00000166428.12" "Expr.ENSG00000149823.8"
## [79] "Expr.ENSG00000101224.17" "Expr.ENSG00000103264.17"
## [81] "Expr.ENSG00000181035.13" "Expr.ENSG00000021300.13"
## [83] "Expr.ENSG00000136997.17" "Expr.ENSG00000103202.12"
## [85] "Expr.ENSG00000106123.11" "Expr.ENSG00000020633.18"
## [87] "Expr.ENSG00000083454.21" "Expr.ENSG00000136068.14"
## [89] "Expr.ENSG00000121966.6" "Expr.ENSG00000077454.15"
## [91] "Expr.ENSG00000160072.19" "Expr.ENSG00000106211.8"
## [93] "Expr.ENSG00000130731.15" "Expr.ENSG00000105643.9"
## [95] "Expr.ENSG00000142546.13" "Expr.ENSG00000120915.13"
## [97] "Expr.ENSG00000125731.12" "Expr.ENSG00000125534.9"
## [99] "Expr.ENSG00000167508.11" "Expr.ENSG00000198355.4"
## [101] "Expr.ENSG00000072310.16" "Expr.ENSG00000133134.11"
## [103] "Expr.ENSG00000197540.7" "Expr.ENSG00000167775.10"
## [105] "Expr.ENSG00000110944.8" "Expr.ENSG00000105486.13"
## [107] "Expr.ENSG00000203485.12" "Expr.ENSG00000197530.12"
## [109] "Expr.ENSG00000157514.16" "Expr.ENSG00000059122.16"
## [111] "Expr.ENSG00000167543.15" "Expr.ENSG00000172534.13"
## [113] "Expr.ENSG00000141524.15" "Expr.ENSG00000100151.15"
## [115] "Expr.ENSG00000090006.17" "Expr.ENSG00000176845.12"
## [117] "Expr.ENSG00000005379.15" "Expr.ENSG00000173762.7"
## [119] "Expr.ENSG00000132635.16" "Expr.ENSG00000143167.11"
## [121] "Expr.ENSG00000180902.17" "Expr.ENSG00000100425.18"
## [123] "Expr.ENSG00000128016.5" "Expr.ENSG00000126351.12"
## [125] "Expr.ENSG00000114812.12" "Expr.ENSG00000013725.14"
## [127] "Expr.ENSG00000162714.12" "Expr.ENSG00000161381.13"
## [129] "Expr.ENSG00000160445.10" "Expr.ENSG00000173846.12"
## [131] "Expr.ENSG00000115085.13" "Expr.ENSG00000173272.15"
## [133] "Expr.ENSG00000162729.13" "Expr.ENSG00000172009.14"
## [135] "Expr.ENSG00000141994.15" "Expr.ENSG00000131634.13"
## [137] "Expr.ENSG00000172543.7" "Expr.ENSG00000152082.13"
## [139] "Expr.ENSG00000144655.14" "Expr.ENSG00000100906.10"
## [141] "Expr.ENSG00000070423.17" "Expr.ENSG00000130787.13"
## [143] "Expr.ENSG00000175602.3" "Expr.ENSG00000125520.13"
## [145] "Expr.ENSG00000167895.14" "Expr.ENSG00000184967.6"
## [147] "Expr.ENSG00000081059.19" "Expr.ENSG00000213366.12"
## [149] "Expr.ENSG00000171223.5" "Expr.ENSG00000172005.10"
## [151] "Expr.ENSG00000139190.16" "Expr.ENSG00000213015.8"
## [153] "Expr.ENSG00000110448.10" "Expr.ENSG00000006015.17"
## [155] "Expr.ENSG00000107872.12" "Expr.ENSG00000099958.14"
## [157] "Expr.ENSG00000100100.12" "Expr.ENSG00000166313.18"
## [159] "Expr.ENSG00000132386.10" "Expr.ENSG00000049089.13"
## [161] "Expr.ENSG00000204564.11" "Expr.ENSG00000156804.7"
## [163] "Expr.ENSG00000143772.9" "Expr.ENSG00000160226.15"
## [165] "Expr.ENSG00000110046.12" "Expr.ENSG00000135960.9"
## [167] "Expr.ENSG00000130522.5" "Expr.ENSG00000107742.12"
## [169] "Expr.ENSG00000140691.16" "Expr.ENSG00000090924.14"
## [171] "Expr.ENSG00000076604.14" "Expr.ENSG00000127528.5"
## [173] "Expr.ENSG00000164880.15" "Expr.ENSG00000167664.8"
## [175] "Expr.ENSG00000123689.5" "Expr.ENSG00000105516.10"
## [177] "Expr.ENSG00000197457.9" "Expr.ENSG00000185022.11"
## [179] "Expr.ENSG00000132819.16" "Expr.ENSG00000034053.14"
## [181] "Expr.ENSG00000168056.15" "Expr.ENSG00000105699.16"
## [183] "Expr.ENSG00000103257.8" "Expr.ENSG00000158106.13"
## [185] "Expr.ENSG00000168209.4" "Expr.ENSG00000174080.10"
## [187] "Expr.ENSG00000109321.10" "Expr.ENSG00000164512.17"
## [189] "Expr.ENSG00000167106.11" "Expr.ENSG00000126368.5"
## [191] "Expr.ENSG00000179094.15" "Expr.ENSG00000158050.4"
## [193] "Expr.ENSG00000163874.10" "Expr.ENSG00000196781.14"
## [195] "Expr.ENSG00000182809.10" "Expr.ENSG00000161835.10"
## [197] "Expr.ENSG00000105655.18" "Expr.ENSG00000177595.17"
## [199] "Expr.ENSG00000118503.14" "Expr.ENSG00000178038.16"
b= weights$ab
ggplot(b, aes(x=Factor5)) +
geom_density()

c= weights$cytokineL
ggplot(c, aes(x=Factor5)) +
geom_density()

d= weights$cell_freq
ggplot(d, aes(x=Factor5)) +
geom_density()

enrichment
against Hallmark
MOFAobject3 = MOFAobject
ens_map = read_tsv('~/Documents/Projects/CMI-PB2/CMI-PB-Oct2023-FINAL/data/raw_prediction_dataset/gene_90_38_export.tsv')
## Rows: 63971 Columns: 12
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (6): biotype, versioned_ensembl_gene_id, description, display_label, syn...
## dbl (6): seq_region_strand, seq_region_start, seq_region_end, canonical_tran...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(ens_map)
## # A tibble: 6 × 12
## biotype seq_region_strand seq_region_start seq_region_end
## <chr> <dbl> <dbl> <dbl>
## 1 miRNA -1 55372940 55373034
## 2 snRNA -1 59450673 59450772
## 3 miRNA 1 41691585 41691678
## 4 protein_coding 1 54201473 54208260
## 5 misc_RNA 1 10287413 10287482
## 6 snRNA 1 109992325 109992431
## # ℹ 8 more variables: versioned_ensembl_gene_id <chr>, description <chr>,
## # canonical_transcript_id <dbl>, display_label <chr>, synonym <chr>,
## # dbprimary_acc <dbl>, gene_id <dbl>, gene_summary <chr>
genes = as.character(rownames(MOFAobject3@data$geneExp$group1))
genes = gsub('Expr.', '', genes)
table(genes %in% ens_map$versioned_ensembl_gene_id)
##
## TRUE
## 6650
ens_map = ens_map[!duplicated(ens_map$versioned_ensembl_gene_id), ]
rownames(ens_map) = ens_map$versioned_ensembl_gene_id
## Warning: Setting row names on a tibble is deprecated.
genes_sym = ens_map[genes, 'display_label']$display_label
rownames(MOFAobject3@data$geneExp$group1) = genes_sym
names(MOFAobject3@intercepts$geneExp$group1) = genes_sym
MOFAobject3@features_metadata$feature[1:6650] = genes_sym
rownames(MOFAobject3@expectations$W$geneExp) = genes_sym
hallmark = fgsea::gmtPathways('~/Documents/References/Genesets/MSigDB/h.all.v2023.2.Hs.symbols.gmt.txt')
# Get all unique genes
all_genes <- sort(unique(unlist(hallmark)))
# Create the binary matrix
binary_matrix <- sapply(all_genes, function(g) {
as.numeric(sapply(hallmark, function(x) g %in% x))
})
# Transpose and convert to a data frame for better visualization
#binary_matrix <- as.data.frame(t(binary_matrix))
binary_matrix[1:5, 1:5]
## A2M AAAS AADAT AARS1 ABAT
## [1,] 0 0 0 0 0
## [2,] 0 0 0 1 0
## [3,] 0 0 0 0 0
## [4,] 0 0 0 0 0
## [5,] 0 0 0 0 0
# Add row names
rownames(binary_matrix) <- names(hallmark)
# Display the binary matrix
binary_matrix[1:5, 1:10]
## A2M AAAS AADAT AARS1 ABAT ABCA1 ABCA2 ABCA3 ABCA4
## HALLMARK_ADIPOGENESIS 0 0 0 0 0 1 0 0 0
## HALLMARK_ALLOGRAFT_REJECTION 0 0 0 1 0 0 0 0 0
## HALLMARK_ANDROGEN_RESPONSE 0 0 0 0 0 0 0 0 0
## HALLMARK_ANGIOGENESIS 0 0 0 0 0 0 0 0 0
## HALLMARK_APICAL_JUNCTION 0 0 0 0 0 0 0 0 0
## ABCA5
## HALLMARK_ADIPOGENESIS 0
## HALLMARK_ALLOGRAFT_REJECTION 0
## HALLMARK_ANDROGEN_RESPONSE 0
## HALLMARK_ANGIOGENESIS 0
## HALLMARK_APICAL_JUNCTION 0
dim(binary_matrix)
## [1] 50 4384
enrichment3.parametric.pos <- run_enrichment(MOFAobject3,
view = "geneExp", factors = 5,
feature.sets = as.matrix(binary_matrix),
sign = "positive",
statistical.test = "parametric"
)
## Intersecting features names in the model and the gene set annotation results in a total of 2153 features.
##
## Running feature set Enrichment Analysis with the following options...
## View: geneExp
## Number of feature sets: 48
## Set statistic: mean.diff
## Statistical test: parametric
## Subsetting weights with positive sign
##
plot_enrichment(enrichment3.parametric.pos,
factor = 'Factor5',
max.pathways = 15
)

plot_enrichment_detailed(enrichment3.parametric.pos,
factor = 'Factor5',
max.genes = 8,
max.pathways = 5
)
## Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

enrichment3.parametric.neg <- run_enrichment(MOFAobject3,
view = "geneExp", factors = 5,
feature.sets = as.matrix(binary_matrix),
sign = "negative",
statistical.test = "parametric"
)
## Intersecting features names in the model and the gene set annotation results in a total of 2153 features.
##
## Running feature set Enrichment Analysis with the following options...
## View: geneExp
## Number of feature sets: 48
## Set statistic: mean.diff
## Statistical test: parametric
## Subsetting weights with negative sign
##
against BTM
btm = fgsea::gmtPathways('~/Documents/References/Genesets/BTM/BTM_for_GSEA_20131008.gmt')
lapply(btm, length)
## $`targets of FOSL1/2 (M0)`
## [1] 12
##
## $`integrin cell surface interactions (I) (M1.0)`
## [1] 29
##
## $`integrin cell surface interactions (II) (M1.1)`
## [1] 12
##
## $`extracellular matrix (I) (M2.0)`
## [1] 30
##
## $`extracellular matrix (II) (M2.1)`
## [1] 45
##
## $`extracellular matrix (III) (M2.2)`
## [1] 14
##
## $`regulation of signal transduction (M3)`
## [1] 47
##
## $`cell cycle and transcription (M4.0)`
## [1] 335
##
## $`cell cycle (I) (M4.1)`
## [1] 145
##
## $`PLK1 signaling events (M4.2)`
## [1] 34
##
## $`myeloid cell enriched receptors and transporters (M4.3)`
## [1] 31
##
## $`mitotic cell cycle - DNA replication (M4.4)`
## [1] 30
##
## $`mitotic cell cycle in stimulated CD4 T cells (M4.5)`
## [1] 35
##
## $`cell division in stimulated CD4 T cells (M4.6)`
## [1] 20
##
## $`mitotic cell cycle (M4.7)`
## [1] 21
##
## $`cell division - E2F transcription network (M4.8)`
## [1] 19
##
## $`mitotic cell cycle in stimulated CD4 T cells (M4.9)`
## [1] 16
##
## $`cell cycle (II) (M4.10)`
## [1] 14
##
## $`mitotic cell cycle in stimulated CD4 T cells (M4.11)`
## [1] 12
##
## $`C-MYC transcriptional network (M4.12)`
## [1] 12
##
## $`cell junction (GO) (M4.13)`
## [1] 11
##
## $`Rho GTPase cycle (M4.14)`
## [1] 10
##
## $`enriched in monocytes (I) (M4.15)`
## [1] 11
##
## $`regulation of antigen presentation and immune response (M5.0)`
## [1] 81
##
## $`T cell activation and signaling (M5.1)`
## [1] 25
##
## $`mitotic cell division (M6)`
## [1] 32
##
## $`enriched in T cells (I) (M7.0)`
## [1] 62
##
## $`T cell activation (I) (M7.1)`
## [1] 50
##
## $`enriched in NK cells (I) (M7.2)`
## [1] 49
##
## $`T cell activation (II) (M7.3)`
## [1] 31
##
## $`T cell activation (III) (M7.4)`
## [1] 14
##
## $`E2F transcription factor network (M8)`
## [1] 18
##
## $`B cell development (M9)`
## [1] 11
##
## $`E2F1 targets (Q3) (M10.0)`
## [1] 33
##
## $`E2F1 targets (Q4) (M10.1)`
## [1] 21
##
## $`enriched in monocytes (II) (M11.0)`
## [1] 189
##
## $`blood coagulation (M11.1)`
## [1] 22
##
## $`formyl peptide receptor mediated neutrophil response (M11.2)`
## [1] 10
##
## $`CD28 costimulation (M12)`
## [1] 10
##
## $`innate activation by cytosolic DNA sensing (M13)`
## [1] 11
##
## $`T cell differentiation (M14)`
## [1] 12
##
## $`Ran mediated mitosis (M15)`
## [1] 13
##
## $`TLR and inflammatory signaling (M16)`
## [1] 45
##
## $`Hox cluster I (M17.0)`
## [1] 16
##
## $`Hox cluster II (M17.1)`
## [1] 12
##
## $`Hox cluster III (M17.2)`
## [1] 12
##
## $`Hox cluster IV (M17.3)`
## [1] 10
##
## $`T cell differentiation via ITK and PKC (M18)`
## [1] 11
##
## $`T cell differentiation (Th2) (M19)`
## [1] 17
##
## $`AP-1 transcription factor network (M20)`
## [1] 15
##
## $`cell adhesion (lymphocyte homing) (M21)`
## [1] 10
##
## $`mismatch repair (I) (M22.0)`
## [1] 27
##
## $`mismatch repair (II) (M22.1)`
## [1] 13
##
## $`RA, WNT, CSF receptors network (monocyte) (M23)`
## [1] 12
##
## $`cell activation (IL15, IL23, TNF) (M24)`
## [1] 15
##
## $`TLR8-BAFF network (M25)`
## [1] 10
##
## $`TBA (M26.0)`
## [1] 30
##
## $`TBA (M26.1)`
## [1] 22
##
## $`TBA (M26.2)`
## [1] 14
##
## $`chemokine cluster (I) (M27.0)`
## [1] 26
##
## $`chemokine cluster (II) (M27.1)`
## [1] 15
##
## $`antigen presentation (lipids and proteins) (M28)`
## [1] 11
##
## $`proinflammatory cytokines and chemokines (M29)`
## [1] 10
##
## $`cell movement, Adhesion & Platelet activation (M30)`
## [1] 20
##
## $`cell cycle and growth arrest (M31)`
## [1] 12
##
## $`platelet activation (I) (M32.0)`
## [1] 23
##
## $`platelet activation (II) (M32.1)`
## [1] 21
##
## $`CORO1A-DEF6 network (I) (M32.2)`
## [1] 20
##
## $`KLF12 targets network (M32.3)`
## [1] 17
##
## $`CORO1A-DEF6 network (II) (M32.4)`
## [1] 15
##
## $`TBA (M32.5)`
## [1] 15
##
## $`TBA (M32.6)`
## [1] 13
##
## $`TBA (M32.7)`
## [1] 11
##
## $`cytoskeletal remodeling (M32.8)`
## [1] 10
##
## $`inflammatory response (M33)`
## [1] 11
##
## $`cytoskeletal remodeling (enriched for SRF targets) (M34)`
## [1] 10
##
## $`signaling in T cells (I) (M35.0)`
## [1] 14
##
## $`signaling in T cells (II) (M35.1)`
## [1] 11
##
## $`T cell surface, activation (M36)`
## [1] 12
##
## $`immune activation - generic cluster (M37.0)`
## [1] 347
##
## $`enriched in neutrophils (I) (M37.1)`
## [1] 52
##
## $`endoplasmic reticulum (M37.2)`
## [1] 19
##
## $`cell division (M37.3)`
## [1] 10
##
## $`chemokines and receptors (M38)`
## [1] 11
##
## $`integrin mediated leukocyte migration (M39)`
## [1] 10
##
## $`complement and other receptors in DCs (M40)`
## [1] 13
##
## $`TBA (M41.0)`
## [1] 17
##
## $`TBA (M41.1)`
## [1] 15
##
## $`TBA (M41.2)`
## [1] 14
##
## $`TBA (M41.3)`
## [1] 10
##
## $`ATF targets network (M41.4)`
## [1] 10
##
## $`platelet activation (III) (M42)`
## [1] 10
##
## $`myeloid, dendritic cell activation via NFkB (I) (M43.0)`
## [1] 15
##
## $`myeloid, dendritic cell activation via NFkB (II) (M43.1)`
## [1] 14
##
## $`T cell signaling and costimulation (M44)`
## [1] 10
##
## $`leukocyte activation and migration (M45)`
## [1] 11
##
## $`cell division (stimulated CD4+ T cells) (M46)`
## [1] 28
##
## $`enriched in B cells (I) (M47.0)`
## [1] 53
##
## $`enriched in B cells (II) (M47.1)`
## [1] 43
##
## $`enriched in B cells (III) (M47.2)`
## [1] 37
##
## $`enriched in B cells (IV) (M47.3)`
## [1] 16
##
## $`enriched in B cells (V) (M47.4)`
## [1] 11
##
## $`TBA (M48)`
## [1] 13
##
## $`transcription regulation in cell development (M49)`
## [1] 47
##
## $`CD1 and other DC receptors (M50)`
## [1] 10
##
## $`cell adhesion (M51)`
## [1] 38
##
## $`T cell activation (IV) (M52)`
## [1] 13
##
## $`inflammasome receptors and signaling (M53)`
## [1] 12
##
## $`BCR signaling (M54)`
## [1] 12
##
## $`TBA (M55)`
## [1] 12
##
## $`suppression of MAPK signaling (M56)`
## [1] 12
##
## $`immuregulation - monocytes, T and B cells (M57)`
## [1] 13
##
## $`B cell development/activation (M58)`
## [1] 11
##
## $`CCR1, 7 and cell signaling (M59)`
## [1] 10
##
## $`lymphocyte generic cluster (M60)`
## [1] 17
##
## $`enriched in NK cells (II) (M61.0)`
## [1] 24
##
## $`enriched in NK cells (KIR cluster) (M61.1)`
## [1] 13
##
## $`enriched in NK cells (receptor activation) (M61.2)`
## [1] 16
##
## $`T & B cell development, activation (M62.0)`
## [1] 44
##
## $`enriched for unknown TF motif CTCNANGTGNY (M62.1)`
## [1] 12
##
## $`regulation of localization (GO) (M63)`
## [1] 12
##
## $`enriched in activated dendritic cells/monocytes (M64)`
## [1] 17
##
## $`IL2, IL7, TCR network (M65)`
## [1] 12
##
## $`TBA (M66)`
## [1] 17
##
## $`activated dendritic cells (M67)`
## [1] 12
##
## $`RIG-1 like receptor signaling (M68)`
## [1] 10
##
## $`enriched in B cells (VI) (M69)`
## [1] 20
##
## $`TBA (M70.0)`
## [1] 20
##
## $`TBA (M70.1)`
## [1] 10
##
## $`enriched in antigen presentation (I) (M71)`
## [1] 18
##
## $`TBA (M72.0)`
## [1] 24
##
## $`TBA (M72.1)`
## [1] 17
##
## $`TBA (M72.2)`
## [1] 12
##
## $`enriched in monocytes (III) (M73)`
## [1] 12
##
## $`transcriptional targets of glucocorticoid receptor (M74)`
## [1] 12
##
## $`antiviral IFN signature (M75)`
## [1] 22
##
## $`DNA repair (M76)`
## [1] 22
##
## $`collagen, TGFB family et al (M77)`
## [1] 35
##
## $`myeloid cell cytokines, metallopeptidases and laminins (M78)`
## [1] 11
##
## $`TBA (M79)`
## [1] 10
##
## $`TBA (M80)`
## [1] 20
##
## $`enriched in myeloid cells and monocytes (M81)`
## [1] 35
##
## $`signal transduction, plasma membrane (M82)`
## [1] 13
##
## $`enriched in naive and memory B cells (M83)`
## [1] 10
##
## $`integrins and cell adhesion (M84)`
## [1] 10
##
## $`platelet activation and degranulation (M85)`
## [1] 12
##
## $`chemokines and inflammatory molecules in myeloid cells (M86.0)`
## [1] 19
##
## $`proinflammatory dendritic cell, myeloid cell response (M86.1)`
## [1] 13
##
## $`transmembrane transport (I) (M87)`
## [1] 24
##
## $`leukocyte migration (M88.0)`
## [1] 51
##
## $`enriched in hepatocyte nuclear factors (I) (M88.1)`
## [1] 13
##
## $`enriched in hepatocyte nuclear factors (II) (M88.2)`
## [1] 12
##
## $`putative targets of PAX3 (M89.0)`
## [1] 16
##
## $`putative targets of PAX3 (M89.1)`
## [1] 10
##
## $`TBA (M90)`
## [1] 12
##
## $`adhesion and migration, chemotaxis (M91)`
## [1] 13
##
## $`lipid metabolism, endoplasmic reticulum (M92)`
## [1] 10
##
## $`TBA (M93)`
## [1] 10
##
## $`growth factor induced, enriched in nuclear receptor subfamily 4 (M94)`
## [1] 12
##
## $`enriched in antigen presentation (II) (M95.0)`
## [1] 24
##
## $`enriched in antigen presentation (III) (M95.1)`
## [1] 15
##
## $`Hox cluster V (M96)`
## [1] 10
##
## $`enriched for SMAD2/3 signaling (M97)`
## [1] 10
##
## $`TBA (M98.0)`
## [1] 18
##
## $`TBA (M98.1)`
## [1] 10
##
## $`TBA (M99)`
## [1] 10
##
## $`MAPK, RAS signaling (M100)`
## [1] 10
##
## $`phosphatidylinositol signaling system (M101)`
## [1] 14
##
## $`TBA (M102)`
## [1] 12
##
## $`cell cycle (III) (M103)`
## [1] 51
##
## $`TBA (M104)`
## [1] 10
##
## $`TBA (M105)`
## [1] 18
##
## $`nuclear pore complex (M106.0)`
## [1] 17
##
## $`nuclear pore complex (mitosis) (M106.1)`
## [1] 11
##
## $`Hox cluster VI (M107)`
## [1] 11
##
## $`TBA (M108)`
## [1] 11
##
## $`receptors, cell migration (M109)`
## [1] 15
##
## $`axon guidance (M110)`
## [1] 10
##
## $`viral sensing & immunity; IRF2 targets network (I) (M111.0)`
## [1] 17
##
## $`viral sensing & immunity; IRF2 targets network (II) (M111.1)`
## [1] 11
##
## $`complement activation (I) (M112.0)`
## [1] 17
##
## $`complement activation (II) (M112.1)`
## [1] 10
##
## $`golgi membrane (I) (M113)`
## [1] 10
##
## $`TBA (M114.0)`
## [1] 39
##
## $`glycerophospholipid metabolism (M114.1)`
## [1] 11
##
## $`cytokines - recepters cluster (M115)`
## [1] 11
##
## $`TBA (M116)`
## [1] 13
##
## $`cell adhesion (GO) (M117)`
## [1] 21
##
## $`enriched in monocytes (IV) (M118.0)`
## [1] 57
##
## $`enriched in monocytes (surface) (M118.1)`
## [1] 17
##
## $`enriched in activated dendritic cells (I) (M119)`
## [1] 11
##
## $`TBA (M120)`
## [1] 11
##
## $`TBA (M121)`
## [1] 12
##
## $`enriched for cell migration (M122)`
## [1] 11
##
## $`enriched in B cell differentiation (M123)`
## [1] 13
##
## $`enriched in membrane proteins (M124)`
## [1] 19
##
## $`TBA (M125)`
## [1] 11
##
## $`double positive thymocytes (M126)`
## [1] 13
##
## $`type I interferon response (M127)`
## [1] 12
##
## $`TBA (M128)`
## [1] 12
##
## $`inositol phosphate metabolism (M129)`
## [1] 13
##
## $`enriched in G-protein coupled receptors (M130)`
## [1] 10
##
## $`TBA (M131)`
## [1] 14
##
## $`recruitment of neutrophils (M132)`
## [1] 11
##
## $`cell adhesion, membrane (M133.0)`
## [1] 14
##
## $`cell cell adhesion (M133.1)`
## [1] 11
##
## $`Membrane, ER proteins (M134)`
## [1] 17
##
## $`enriched in plasma membrane proteins (I) (M135.0)`
## [1] 16
##
## $`enriched in plasma membrane proteins (II) (M135.1)`
## [1] 15
##
## $`TBA (M136)`
## [1] 17
##
## $`TBA (M137)`
## [1] 16
##
## $`enriched for ubiquitination (M138)`
## [1] 11
##
## $`lysosomal/endosomal proteins (M139)`
## [1] 11
##
## $`extracellular matrix, complement (M140)`
## [1] 14
##
## $`TBA (M141)`
## [1] 27
##
## $`transmembrane and ion transporters (I) (M142)`
## [1] 13
##
## $`nuclear pore, transport; mRNA splicing, processing (M143)`
## [1] 11
##
## $`cell cycle, ATP binding (M144)`
## [1] 16
##
## $`cytoskeleton/actin (SRF transcription targets) (M145.0)`
## [1] 16
##
## $`cytoskeleton/actin (SRF transcription targets) (M145.1)`
## [1] 14
##
## $`MHC-TLR7-TLR8 cluster (M146)`
## [1] 17
##
## $`intracellular transport (M147)`
## [1] 17
##
## $`TBA (M148)`
## [1] 11
##
## $`TBA (M149)`
## [1] 10
##
## $`innate antiviral response (M150)`
## [1] 12
##
## $`TBA (M151)`
## [1] 10
##
## $`TBA (source: B cells) (M152.0)`
## [1] 25
##
## $`TBA (source: naive B cells) (M152.1)`
## [1] 21
##
## $`TBA (source: memory B cells) (M152.2)`
## [1] 18
##
## $`TBA (M153)`
## [1] 16
##
## $`amino acid metabolishm and transport (M154.0)`
## [1] 33
##
## $`transmembrane transport (SLC cluster) (M154.1)`
## [1] 17
##
## $`G protein coupled receptors cluster (M155)`
## [1] 10
##
## $`plasma cells & B cells, immunoglobulins (M156.0)`
## [1] 43
##
## $`plasma cells, immunoglobulins (M156.1)`
## [1] 36
##
## $`enriched in NK cells (III) (M157)`
## [1] 14
##
## $`interferon alpha response (I) (M158.0)`
## [1] 16
##
## $`interferon alpha response (II) (M158.1)`
## [1] 13
##
## $`G protein mediated calcium signaling (M159)`
## [1] 10
##
## $`leukocyte differentiation (M160)`
## [1] 16
##
## $`TBA (M161)`
## [1] 14
##
## $`plasma membrane, cell junction (M162.0)`
## [1] 18
##
## $`cell junction (M162.1)`
## [1] 11
##
## $`enriched in neutrophils (II) (M163)`
## [1] 14
##
## $`xenobiotic metabolism (M164)`
## [1] 10
##
## $`enriched in activated dendritic cells (II) (M165)`
## [1] 35
##
## $`TBA (M166)`
## [1] 17
##
## $`enriched in cell cycle (M167)`
## [1] 15
##
## $`enriched in dendritic cells (M168)`
## [1] 19
##
## $`mitosis (TF motif CCAATNNSNNNGCG) (M169)`
## [1] 16
##
## $`TBA (M170)`
## [1] 12
##
## $`heme biosynthesis (I) (M171)`
## [1] 11
##
## $`enriched for TF motif TTCNRGNNNNTTC (M172)`
## [1] 10
##
## $`erythrocyte differentiation (M173)`
## [1] 14
##
## $`TBA (M174)`
## [1] 24
##
## $`cell development (M175)`
## [1] 13
##
## $`TBA (M176)`
## [1] 10
##
## $`TBA (M177.0)`
## [1] 34
##
## $`TBA (M177.1)`
## [1] 15
##
## $`enriched for promoter motif NATCACGTGAY (putative SREBF1 targets) (M178)`
## [1] 10
##
## $`enriched for TF motif PAX3 (M179)`
## [1] 10
##
## $`TBA (M180)`
## [1] 11
##
## $`nucleotide metabolism (M181)`
## [1] 10
##
## $`enriched in DNA interacting proteins (M182)`
## [1] 16
##
## $`TBA (M183)`
## [1] 10
##
## $`TBA (M184.0)`
## [1] 18
##
## $`TBA (M184.1)`
## [1] 11
##
## $`TBA (M185)`
## [1] 14
##
## $`TBA (M186)`
## [1] 11
##
## $`metabolism in mitochondria, peroxisome (M187)`
## [1] 12
##
## $`TBA (M188)`
## [1] 10
##
## $`extracellular region cluster (GO) (M189)`
## [1] 16
##
## $`TBA (M190)`
## [1] 11
##
## $`transmembrane transport (II) (M191)`
## [1] 18
##
## $`TBA (M192)`
## [1] 12
##
## $`TBA (M193)`
## [1] 11
##
## $`TBA (M194)`
## [1] 14
##
## $`muscle contraction, SRF targets (M195)`
## [1] 12
##
## $`platelet activation - actin binding (M196)`
## [1] 17
##
## $`TBA (M197)`
## [1] 14
##
## $`TBA (M198)`
## [1] 12
##
## $`platelet activation & blood coagulation (M199)`
## [1] 13
##
## $`antigen processing and presentation (M200)`
## [1] 11
##
## $`TBA (M201)`
## [1] 18
##
## $`enriched in extracellular matrix & associated proteins (M202)`
## [1] 12
##
## $`TBA (M203)`
## [1] 13
##
## $`chaperonin mediated protein folding (I) (M204.0)`
## [1] 14
##
## $`chaperonin mediated protein folding (II) (M204.1)`
## [1] 10
##
## $`TBA (M205)`
## [1] 11
##
## $`Wnt signaling pathway (M206)`
## [1] 14
##
## $`TBA (M207)`
## [1] 10
##
## $`TBA (M208)`
## [1] 10
##
## $`lysosome (M209)`
## [1] 10
##
## $`extracellular matrix, collagen (M210)`
## [1] 32
##
## $`TBA (M211)`
## [1] 11
##
## $`purine nucleotide biosynthesis (M212)`
## [1] 12
##
## $`regulation of transcription, transcription factors (M213)`
## [1] 21
##
## $`TBA (M214)`
## [1] 12
##
## $`small GTPase mediated signal transduction (M215)`
## [1] 16
##
## $`respiratory electron transport chain (mitochondrion) (M216)`
## [1] 12
##
## $`TBA (source: B cells) (M217)`
## [1] 10
##
## $`TBA (M218)`
## [1] 15
##
## $`respiratory electron transport chain (mitochondrion) (M219)`
## [1] 18
##
## $`TBA (M220)`
## [1] 10
##
## $`TBA (M221)`
## [1] 16
##
## $`heme biosynthesis (II) (M222)`
## [1] 13
##
## $`enriched in T cells (II) (M223)`
## [1] 13
##
## $`transmembrane and ion transporters (II) (M224)`
## [1] 11
##
## $`metabolism of steroids (M225)`
## [1] 10
##
## $`proteasome (M226)`
## [1] 12
##
## $`translation initiation (M227)`
## [1] 10
##
## $`olfactory receptors (M228)`
## [1] 13
##
## $`TBA (M229)`
## [1] 11
##
## $`cell cycle, mitotic phase (M230)`
## [1] 12
##
## $`respiratory electron transport chain (mitochondrion) (M231)`
## [1] 11
##
## $`enriched for TF motif TNCATNTCCYR (M232)`
## [1] 13
##
## $`TBA (M233)`
## [1] 15
##
## $`transcription elongation, RNA polymerase II (M234)`
## [1] 14
##
## $`mitochondrial cluster (M235)`
## [1] 19
##
## $`TBA (M236)`
## [1] 14
##
## $`golgi membrane (II) (M237)`
## [1] 17
##
## $`respiratory electron transport chain (mitochondrion) (M238)`
## [1] 17
##
## $`enriched in calcium signaling proteins (M239)`
## [1] 15
##
## $`chromosome Y linked (M240)`
## [1] 17
##
## $`TBA (M241)`
## [1] 10
##
## $`TBA (M242)`
## [1] 11
##
## $`TBA (M243)`
## [1] 12
##
## $`TBA (M244)`
## [1] 13
##
## $`translation initiation factor 3 complex (M245)`
## [1] 13
##
## $`TBA (M246)`
## [1] 16
##
## $`enriched in nuclear pore complex interacting proteins (M247)`
## [1] 14
##
## $`TBA (M248)`
## [1] 11
##
## $`TBA (M249)`
## [1] 10
##
## $`spliceosome (M250)`
## [1] 12
##
## $`T cell surface signature (S0)`
## [1] 27
##
## $`NK cell surface signature (S1)`
## [1] 48
##
## $`B cell surface signature (S2)`
## [1] 169
##
## $`Plasma cell surface signature (S3)`
## [1] 24
##
## $`Monocyte surface signature (S4)`
## [1] 94
##
## $`DC surface signature (S5)`
## [1] 82
##
## $`CD4 T cell surface signature Th1-stimulated (S6)`
## [1] 9
##
## $`CD4 T cell surface signature Th2-stimulated (S7)`
## [1] 13
##
## $`Naive B cell surface signature (S8)`
## [1] 54
##
## $`Memory B cell surface signature (S9)`
## [1] 37
##
## $`Resting dendritic cell surface signature (S10)`
## [1] 75
##
## $`Activated (LPS) dendritic cell surface signature (S11)`
## [1] 37
# Get all unique genes
all_genes <- sort(unique(unlist(btm)))
# Create the binary matrix
binary_matrix_btm <- sapply(all_genes, function(g) {
as.numeric(sapply(btm, function(x) g %in% x))
})
# Transpose and convert to a data frame for better visualization
#binary_matrix_btm <- as.data.frame(t(binary_matrix_btm))
binary_matrix_btm[1:5, 1:5]
## A1CF A2BP1 A2M ABCA1 ABCA13
## [1,] 0 0 0 0 0
## [2,] 0 0 0 0 0
## [3,] 0 0 0 0 0
## [4,] 0 0 0 0 0
## [5,] 0 0 0 0 0
# Add row names
rownames(binary_matrix_btm) <- names(btm)
# Display the binary matrix
binary_matrix_btm[1:5, 1:10]
## A1CF A2BP1 A2M ABCA1 ABCA13
## targets of FOSL1/2 (M0) 0 0 0 0 0
## integrin cell surface interactions (I) (M1.0) 0 0 0 0 0
## integrin cell surface interactions (II) (M1.1) 0 0 0 0 0
## extracellular matrix (I) (M2.0) 0 0 0 0 0
## extracellular matrix (II) (M2.1) 0 0 0 0 0
## ABCA6 ABCA8 ABCB4 ABCB6 ABCC13
## targets of FOSL1/2 (M0) 0 0 0 0 0
## integrin cell surface interactions (I) (M1.0) 0 0 0 0 0
## integrin cell surface interactions (II) (M1.1) 0 0 0 0 0
## extracellular matrix (I) (M2.0) 0 0 0 0 0
## extracellular matrix (II) (M2.1) 0 0 0 0 0
dim(binary_matrix)
## [1] 50 4384
enrichment4.parametric.pos <- run_enrichment(MOFAobject3,
view = "geneExp", factors = 5,
feature.sets = as.matrix(binary_matrix_btm),
sign = "positive",
statistical.test = "parametric"
)
## Intersecting features names in the model and the gene set annotation results in a total of 1277 features.
##
## Running feature set Enrichment Analysis with the following options...
## View: geneExp
## Number of feature sets: 124
## Set statistic: mean.diff
## Statistical test: parametric
## Subsetting weights with positive sign
##
plot_enrichment(enrichment4.parametric.pos,
factor = 'Factor5',
max.pathways = 15
)

plot_enrichment_detailed(enrichment4.parametric.pos,
factor = 'Factor5',
max.genes = 8,
max.pathways = 5
)
## Warning: ggrepel: 3 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

enrichment4.parametric.neg <- run_enrichment(MOFAobject3,
view = "geneExp", factors = 5,
feature.sets = as.matrix(binary_matrix_btm),
sign = "negative",
statistical.test = "parametric"
)
## Intersecting features names in the model and the gene set annotation results in a total of 1277 features.
##
## Running feature set Enrichment Analysis with the following options...
## View: geneExp
## Number of feature sets: 124
## Set statistic: mean.diff
## Statistical test: parametric
## Subsetting weights with negative sign
##
plot_enrichment(enrichment4.parametric.neg,
factor = 'Factor5',
max.pathways = 15
)

plot_enrichment_detailed(enrichment4.parametric.neg,
factor = 'Factor5',
max.genes = 8,
max.pathways = 5
)
## Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

other modules
enrichment.parametric <- run_enrichment(MOFAobject3,
view = "geneExp", factors = c(1:15),
feature.sets = as.matrix(binary_matrix_btm),
sign = "positive",
statistical.test = "parametric"
)
## Intersecting features names in the model and the gene set annotation results in a total of 1277 features.
##
## Running feature set Enrichment Analysis with the following options...
## View: geneExp
## Number of feature sets: 124
## Set statistic: mean.diff
## Statistical test: parametric
## Subsetting weights with positive sign
##
weights <- get_weights(MOFAobject3, views = "all", factors = "Factor5", scale = F)
weights.scale <- get_weights(MOFAobject3, views = "geneExp", factors = "Factor5", scale = T)
a= weights$ab
ggplot(a, aes(x=Factor5)) +
geom_density()

a['IgG_PT', ]
## [1] -8.075284e-07
#weights.scale$geneExp[c('CCL3', 'CXCL8', 'IL1B', 'CCL4'), ]
a = a[order(a[,1]), ]
head(a)
## IgG1_TT IgG1_DT IgG1_PRN IgG3_FHA IgG1_FHA
## -0.0013979687 -0.0008164969 -0.0005639915 -0.0002660241 -0.0002602516
## IgG3_DT
## -0.0002376422
tail(a)
## IgG2_FHA IgG3_PT IgG4_TT IgG2_OVA IgG4_DT IgG4_OVA
## 0.0001782034 0.0001866793 0.0002234760 0.0003458956 0.0003545961 0.0006920662
d = as.data.frame(tail(names(a), n=200))
#write_tsv(d, file='data/CCL3_MOFA_F3.tsv', col_names = T)
d
## tail(names(a), n = 200)
## 1 IgG1_TT
## 2 IgG1_DT
## 3 IgG1_PRN
## 4 IgG3_FHA
## 5 IgG1_FHA
## 6 IgG3_DT
## 7 IgG4_FHA
## 8 IgG_PRN
## 9 IgG2_DT
## 10 IgG4_PRN
## 11 IgG3_PRN
## 12 IgG1_PT
## 13 IgG2_PRN
## 14 IgG3_OVA
## 15 IgG_PT
## 16 IgG4_PT
## 17 IgG3_FIM2/3
## 18 IgG4_FIM2/3
## 19 IgG2_FIM2/3
## 20 IgG_FHA
## 21 IgG1_OVA
## 22 IgG2_PT
## 23 IgG3_TT
## 24 IgG1_FIM2/3
## 25 IgG2_TT
## 26 IgG2_FHA
## 27 IgG3_PT
## 28 IgG4_TT
## 29 IgG2_OVA
## 30 IgG4_DT
## 31 IgG4_OVA
cytokineLs
plot_top_weights(MOFAobject3,
view = "cytokineL",
factor = 5,
nfeatures = 5
)

plot_top_weights(MOFAobject3,
view = "cell_freq",
factor = 5,
nfeatures = 5
)

plot_data_scatter(MOFAobject,
view = "ab", # view of interest
factor = 5, # factor of interest
features = 5, # number of features to plot (they are selected by weight)
add_lm = TRUE, # add linear regression
color_by = "Meta.infancy_vac"
)
